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Abstract. - We propose a scheme to measure the mass of a single particle using the nonlinear 
response of a 2D nanoresonator with degenerate eigenmodes. Using numerical and analytical cal- 
culations, we show that by driving a square graphene nanoresonator into the nonlinear regime, 
simultaneous determination of the mass and position of an added particle is possible. More- 
over, this scheme only requires measurements in a narrow frequency band near the fundamental 
resonance. 



Introduction. — Nanoelectromechanical (NEM) res- 
onators hold promise as ultrasensitive mass detectors p~lf2] . 
NEM mass sensors (NEM-MS) rely on a resonant fre- 
quency shift Aid due to an added mass AM. However, 
as opposed to detecting a single adsorbed particle, to ac- 
tually measure its mass AM from Aw, the position of the 
particle must be known. Proposed position determination 
schemes 3-6 rely on detectors to measure the frequency 
shifts of several vibration modes. While this poses no 
problems in principle, it causes practical difficulties for 
NEM-MS operating in the GHz regime. 

We propose a detection scheme that only requires mea- 
surements in a single narrow band centered at the fun- 
damental mode resonance frequency of a square 2D res- 
onator. Our method uses the nonlinear response of the 
resonator by exploiting the interaction between vibration 
modes to make information about higher modes available 
at the fundamental frequency. We illustrate by showing, 
analytically and numerically, how the nonlinear response 
of micrometer-size graphene resonators [Zl[8] can be used 
for single particle mass measurements with zeptogram pre- 
cision at room temperatures. 

Several other technology tracks are being considered for 
NEM-MS devices. One is downscaling of Si-MEMS 
[TTUl4j where the present state-of-the-art give a minimum 
detectable mass of ~ 10 zg [9]. Another track relies on 
carbon nanotubes (CNTs) [TS] and has already reached 
sub-zg levels [T0lfT6rtT8] . However, after the discovery 
of graphene [19], novel 2D NEMS devices have been ex- 



Suspended graphene 



Added mass 




Gate 



Fig. 1: Possible realisation of a NEM mass spectrometer using 
a suspended square graphene sheet with all edges clamped. Be- 
low the graphene an electrostatic gate for actuation and trans- 
duction is placed symmetrically with respect to the X-axis and 
asymmetrically with respect to the Y - axis. By electrostatic 
actuation of vibration modes, a mass AM located at an arbi- 
trary position Xm = (Xm,Ym) can be determined. 



plored [20-23 , including mass detectors with zg sensitiv- 
ity [7] . Apart from increasing the adsorbtion cross-section, 
2D-NEMS can also have degenerate fiexural modes. As we 
show, this degeneracy makes possible to distinguish single- 
particle from multi-particle adsorption. Graphene also 
represents the ultimate material for 2D-NEMS through 
its combination of large strength and low mass. 
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System. — We consider a square graphene sheet with 
mass M and side length Lq suspended in the XF-plane 
above an actuation gate (See Figure [TJ. The sheet is sim- 
ply clamped at all edges. The gate geometry, which has a 
symmetry line parallel to the Y-axis, is chosen such that 
the fundamental and higher order modes can be excited. 
The transverse deflection u>(X, t) of the membrane is given 
by 

pw + cw- Y d d T t 9 S w ) = P z(X,t). (1) 

$=X,Y 

Here P z is the external pressure on the sheet. This pres- 
sure comes from the electric biasing on the gate elec- 
trode. The exact geometry of the gate, and the exact 
X-dependence of P z need not be known. It suffices that P z 
has the proper symmetry. And, Tx = Ty = To + 7i|Vu>| 2 
are sheet tension components where T is an initial ten- 
sion and T\ m 112 N/m. Equation (JTJ is nonlinear due 
to stretching-induced tension |S]. For a particle with rel- 
ative mass e = AM/M adsorbed at Xm, the density is 
p(X) = po + AM(5(X - X M ), where 5(X) is the 2D delta 
function and po is the density of graphene. 

For future convenience, we begin by rescaling Eq.lfTj) 
into a dimensionless form. We do this by introduc- 
ing the length and time scales ho = LoyTo/Tx and 
to = Lq y/ po /Tq , we write the deflection as u(x, t) = 
w(Lo x i to T )/ho. Equation (fTJ) then becomes 

[l+ e( 5(x-XAf)]w+7M-V 2 ii- Y d 6 (\Vu\ 2 d £ u) = p z (2) 

where 7 = ct Q /po and p z = P z t\ / '(p h ) . 

Linear response. — We consider first small deflec- 
tions where Tx,y ~ To, and the resonator is in the linear 
regime. The eigenmodes are then determined from 

-u; 2 [l + e5(x-x M )]u- V 2 u = 0, xe[0,l] 2 . (3) 

Without adsorbed particles e = 0, the first three 
mode shapes are 0io = 2 sin (irx) sin (Try), feo = 
2 sin (2ttx) sin (ny), $30 — 2 sin (nx) sin (2ny), with eigen- 
frcqucncics u\ = 2-zr 2 and w 2 = w| = 57r 2 . To linear or- 
der in e, adding a mass at xm leads to lo\ — w 2 (l — £0 2 ), 
w 2 = <Aq{1 - eA" 2 ), and uj 3 = uj 30 . Here m = mO (x Af ) 
and A" = [0| + 0i] 1 ^ 2 - To zeroth order in e, 0i = 1O , 

02 = [02020 + 0303o]/A" and 3 = [02030 - 03020 ]/ A". 

These solutions are illustrated in Fig. [2j 

For a two-fold degenerate mode, the frequency of one 
mode is lowered due to particle adsorbtion. The other 
mode will not change frequency since it has a nodal line 
passing through the location Xm- This allows a simple 
test to see if more than one particle has been adsorbed. 
A multi-particle adsorption results in frequency shifts for 
both the initially degenerate modes. 

Nonlinear response. — To study the nonlinear dy- 
namics of the system, we expand the scaled deflection u 
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Fig. 2: Amplitudes for the three lowest flexural eigenmodes as 
functions of drive frequency oj for weak driving. Dashed lines: 
Linear response without added mass. The unperturbed mode 
shapes c/uo, (A20 and 4>zo are indicated on the plaquettes where 
the locations of nodelines antinodes are shown. The modes <j>20 
and <j>zo are degenerate. Solid lines: Linear response in the 
presence of an added mass. The mode functions are 0i, 02 
and c/>3 with shapes indicated on the plaquettes. The blue dots 
show the position of the added mass. 

in Eq. ([1]) in the eigenmodes m (x) of the linear problem 
[Eq.© with e ^ 0] as u(x, r) = £^ =1 u m (r)0 m (x). This 
yields a system of coupled Duffing equations for the mode 
amplitudes u m 

00 

D m (u m + 0Jj n u m ) + ■yii m + Y A mrst u r u s u t = p m - (4) 

rst— 1 

Here D m = l+e0 m (x M ) 2 = l + £0 2 „, A mrst = J dx(V0 m - 
V0 r )(V0 s • V0t), and p m = J dx<j) m p z . As e < 1 we have 
to lowest order in e, D" 1 w 1 - e<j> 2 m w w^/i^io- 

(u m + u? m u m ) + y[1 - e0 2 Jw m 

OC 

+ Y A mrst [l - e<j) 2 m ]u r u s ut = p m [l - e0 2 J. (5) 

rst— 1 

In what follows we will consider the weakly nonlinear 
regime. The cubic nonlinearities in Eq. ([5]) can be then be 
treated using the method of averaging (Krylov-Bogoliubov 
method). In this method, both the damping jii, the driv- 
ing p m , and the terms of order u 3 are of the same order and 
small (see for instance Ref. [28]). Formally, 7 can in this 
method be treated as a small parameter of a perturbation 
expansion. To simplify the analysis, terms of order 0(e~f) 
can then be considered as higher order terms and omitted. 
Further, only drive frequencies close to u>io and W20 = w 30 
are used and equations for the three lowest modes suffice. 
These approximations give 

u\ + "fiii + (w 2 + 5[Au2 + Au 2 ])ui + Au\ = p\ 
u 2 + 7it 2 + (w 2 + h[Au\ + Cm 2 ])m2 + Bui = Pi 
U3 + 7"3 + K 2 + 5{Auj + Cu\])u z + Bui = P3 (6) 

where A = 5tt 4 , B = 161tt 4 /4 + 3tt 4 2 2 /(2A' 4 ) and C 
4l7r 4 /5. The ultimate justification for the approximations 
leading up to Eq. ([6]) are the comparisons of the theoretical 
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Fig. 3: Amplitudes for modes 2 and 3 as functions of drive fre- 
quency uj for a square membrane with an added mass. Solid 
lines: Nonlinear response. Dashed lines: Linear response 
(see Fig.[5|. By driving both modes into the nonlinear regime, 
the parameter r [see Eq.[TT] can be obtained from the frequency 
shifts Au>c2 and Alu C 3. The parameter r defines the nodal line 
of mode 3. Both u) C 2 and u> C 3 are measured by sweeping ui 
downwards. Solid curves were obtained by numerical integra- 
tion of Eq. © with a mass fraction AM/M = 0.08% located 
at (x m ,Vm) = (0.81,0.20) (quality factor Qi = 3000). Dash- 
dotted line: Above the frequencies lj C 2,c3 hysteretic behavior 
can be observed by sweeping cj upwards. 

treatment of the system ((5]) with the numerical simulations 
of the full equations ((H) . 

For the external force of the form p z (x, r) = p(r)g(x) 
where g obeys the symmetry relation <?(x) = g(\x— 0.5 1 , y), 
the source terms can be written as 

pi( t ) = D M T ) 

P 2 (t) = D 2 p(T)cos(iry M ) 
P 3 (t) = D 2 p{t)cos(itxm)- 

Here 

Di = 2 / cfoc, sin(7r:r) sin(7ry)g(x) 



and 



Do 



^Jdx sin(7ra;) sin(27ry)(7(x) 
yj cos 2 ttxm + cos 2 TvyM 



In the expressions for for the source terms p n , the form 
of the driving force, <?(x) is included in the coefficients 
Di 2- We again stress that the exact form of <?(x) is not im- 
portant, and need not be known, as long as it has the sym- 
metry property g(x) = g(\x — 0.5|, y). It is this symmetry 
property which causes the same coeffecient D 2 to appear 
in both the source terms p 2 and p 3 . Hence, any measur- 
able quantity which depends only on the ratio p 2 /pi will 
thus be a function of only the particle position xm- This 
will be used in the mass measurment scheme presented 
below. 

Mass measurement. — To determinine the position 
of the adsorbed mass we will use the parameters r and s 



defined as 

r = cos(7n/ M ) 2 / cos(ttiz;m) 2 (7) 
s = 1 - [cos 2 (ttxm) + cos 2 (7ry M )]- (8) 

The quantity s is related to the frequency shifts in the 
linear response regime through 



1 L0i 



10 ujf - uif 



(9) 



This parameter can thus be determined by applying a 
weak harmonic drive of the form p(r) — cos(wt) and mon- 
itoring the location of resonances. Driving the system 
harder, still with a single frequency, puts it in the non- 
linear regime. However, for a single frequency excitation 
in the weakly non-linear regime, the coupling between the 
equations in ^ can be ignored and the system turns into 
three uncoupled Duffing equations. 



u\ + "fill + u\ui + Au\ = px 
■ uj 2 u 2 + Bui = P2 



u 2 

"3 



JU 2 
■JU 3 



w 2 u 3 + Bui = Ps 



(10) 



Characteristic for a driven Duffing oscillator in the non- 
linear regime is the bistability region in parameter space 
where the system oscillates with either small or large am- 
plitude depending on the initial conditions. This leads to 
the characteristic hysteresis loops seen in figure [3] 

The parameter r can be related to the frequency shifts 
by noting that the ratio of the forces p 2 [r) and p 3 (r) in 
Eq. ((6|) is given by \Jr. As shown in appendix, the edges 
of the hysteresis loops depend on the applied forces as 
K 2 „ - W 2 ) 3 « (9/4) 2 £p 2 (n = 2, 3) so that 



J c1 



u c3 



(11) 



Hence, frequency measurements in the linear and nonlin- 
ear regimes can be used to determine r and s. From r and 
s the position of the adsorbed particle can be deduced (up 
to symmetry of the structure). Knowing the position (in 
terms of r and s) allows calculation of the mass respon- 
sivity IZi of the fundamental mode <f>\ by calculating the 
linear frequency shift 



Ki(xm) 1=3 -2wi 



•r)(l 



(1 



(12) 



which gives the added mass AM = eM = K^MAcj!. 

The result presented here rests on three main equations 
© , (ITU and (|T2|) . To obtain this result we have made two 
crucial assumptions relating to the symmetry of the sys- 
tem; the symmetry leading to mode degeneracy and the 
symmetry of the gate. In any real situation, these sym- 
metries will not be exact and it is relevant to question to 
what extent these symmetries will need to be fulfilled. For 
a complete error-analysis, one must analyze the detailed 
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reasons for lifting the degeneracies. While such a detailed 
analysis is beyond the scope of the present work, some ob- 
servations can be readily made. Firstly, the most crucial 
symmetry is that of the membrane. For the scheme pre- 
sented here to be relevant thus puts constraints on the in- 
trinsic mode splitting Aw 2 3 = o; 30 — o; 20 . The first of these 
constraints is AW23 w 20 — w 2- When this inequality is 
fulfilled, the effect of an adsorbed particle on the nearly 
degeneraty modes is larger than the effect of imperfections 
leading to the intrinsic splitting. A second criterion, which 
is less obvious, is that 

Aw 2 3 <C w 3 - w 2 

This criterion means that mode 3 does not shift apprecia- 
bly when the particle is added. 

Narrowband scheme. — Above, we have demon- 
strated that frequency measurements can be used to deter- 
mine the position and mass of the adsorbed particle. We 
now show that, by exploiting the nonlinearities in the sys- 
tem, this information can be obtained by measuring only 
in a narrow frequency band near the fundamental mode 
frequency uj\. 

Equations © represent a system of three coupled Duff- 
ing oscillators for the modes amplitudes u n [n — 1,2,3]. 
Here, the effective resonant frequency of a mode depends 
not only on the oscillation amplitude of the mode itself 
but also on the amplitudes of other modes so that for in- 
stance uj\ increase by approximately 5A^ 23 (it^) where 
(•) denotes time-average over an oscillation period. This 
allows us to choose to use the fundamental mode to mon- 
itor the amplitudes of modes 2 and 3 as follows: In the 
first step, the system is excited with a single frequency 
signal p(r) = pa cos(ojt) and the frequency oji of the fun- 
damental mode in the linear regime is determined. The 
frequency of this excitation, and detection, is henceforth 
kept fixed at u)\. A second excitation signal pscos(wr) 
is superimposed on the signal at frequency uj\. When the 
amplitude ps is low, the excitation of mode 2 in the linear 
regime for ui — UJ2 can be detected as a reduction of the 
oscillation amplitude of the fundamental mode. This is 
because the effective frequency of the fundamental mode 
is shifted away from lu± due to the excitation of mode 2. 
Finally, when ps is increased, the mode 2 is driven into 
the nonlinear regime and w C 2 can be determined. Simi- 
larly, W3 and w C 3 can be obtained. The effect of the mode 
interaction between the fundamental mode and modes 2 
and 3 are shown in Fig. |H 

At first hand one may object to this scheme by noting 
that when the fundamental mode is strongly excited, it af- 
fects the frequencies w 2 and w c2 . However, since both lo\ 
and w 2 2 shift by the same amount, these shifts cancel out 
(to first order) in the expression for r. The cancellation 
occurs also in the expression for s if the resonant frequen- 
cies cj„o before mass adsorption are determined through 
the same narrowband scheme. 
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Fig. 4: Mode amplitudes obtained by numerical integration 
of the system ([4]) using a gate signal pa cos(wit) + ps cos(tJi). 
Upper panel: Amplitude of mode 1 as function of variable 
drive frequency cj. Lower panel: Amplitudes for modes 2 
and 3 as functions of drive frequency uj. The frequency wi 
is fixed at the resonance of mode 1 while oj is varied. Due 
to nonlinearity the modes couple. This causes the resonant 
frequency of mode 1 to depend on the amplitudes of modes 2 
and 3. It will thus shift away from Wi for finite amplitudes of 
modes 2 and 3. Hence, by measuring the response of mode 1, 
the responses of modes 2 and 3 can be probed by measuring 
only in a narrow frequency band around uii. 

Mesurement sensitivity and range. Wc now 

consider sensitivity and range. In the NEM-MS experi- 
ments reported in the literature [5HT8]. the sensitivity is 
usually taken as the smallest detectable mass. In our case 
this occurs when the particle is adsorbed at the sweet 
spot of the resonator at xj./ = (0.5,0.5). This leads to 
AM m i n = 0.5(Auji/ujio) m i n M . The intrinsic limitation on 
| Aw/a; | comes from thermomechanical noise that deter- 
mines how small resonance shift can be reliably detected. 
If the detector bandwidth Aw is narrower than the reso- 
nance at ui we have |Aw/w| > Q^^O -011 "/ 20 [12]. Here 
DR n is the dynamic range of mode n and Qi the quality 
factor of the fundamental mode. For modes n = 1, 2, 3 we 
find 



DR n = 10 log 



Rn 


(T \ 


ToL 


.Qi 




k B T\ 



where R\ ps 0.6 and i?2 = R3 ~ 0.3. For a device with 
L a = 1 fim, Qi = 3000 and wi/(2tt) = 2 GHz we find 
AM min fa iMQ^lfr 2 - 5 ps 0.5 zg at T = 300 K. At lower 
temperatures the sensitivity improves as T 1 / 2 . 

Thermal fluctuations also influence the determination 
of the frequencies w c2 , C 3. If the system performs low- 
amplitude oscillations with uj close to oj c , thermal fluc- 
tuations can cause transitions to the high-amplitude state 
before u> c is reached. To accurately determine ui c we must 
have W <C w c where W is the rate for transitions to the 
high amplitude state. This rate obeys W oc e - RE T/(k B T) 
where E T = T 2 ig/Ti [27]. As demonstrated in Ref. [24], 
the strong exponential dependence of W on system pa- 
rameters can for NEMS lead to an enhanced sensitivity 
in the measurements of w C 2, C 3 compared to the frequency 
measurements in the linear regimes. 

We now consider the range of masses that can be reliably 
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Fig. 5: (a) Maximal values of e = AM/M due to limita- 
tions of first order perturbation theory. Within each con- 
tour, mass fractions up to e m ax can be determined with a 
5% accuracy, (b) Contours of minimum eQi where Eq. 
is applicable. E.g., in the shaded area Eq. is valid for 

e > I.6/Q1. (c) Determination of randomly deposited masses 
using numerical integration of Eq. Q for a membrane with 
Qi — 3000. The masses were uniformly distributed in the 
range 0.02% < e < 0.35%. Frequencies were determined us- 
ing an accuracy of [Aw/w| ~ 0.5 ■ 10 -4 . The positions of the 
deposited masses are shown by shaded symbols. The open 
symbols were obtained using Eqs. @ and ((TTJ . The size of 
the markers are proportional to e. The dashed lines indicate 
regions where |(e — e cxact ) /e exa ct | is less than 2% or 10%. 

measured with the nonlinear mass determination scheme 
presented above. This must not be confused with the sen- 
sitivity discussed above which only considers the minimum 
detectable mass change. The range includes both upper 
and lower bounds on e = AM/M. The upper bound arises 
from omitting terms of C(e 2 ) and higher in the relation 
Awi = lZie + (D(e 2 ). Fig. [5^ shows contours on a quadrant 
of the unit square corresponding to the membrane. Each 
contour encloses a region where the relative error due to 
omitting terms of 0(e 2 ) is less than 5%. For instance, 
masses with e up to e max = 0.1% can only be determined 
with a relative error less than 5% if they are located in- 
side the e max = 0.1%-contour. The upper bound can be 
improved upon by using numerically calculated values of 
Awi(c,xm) instead of perturbation theory. 

Specific to this scheme is that to determine r in Eq. (ITTT) , 
the regions of multivalued response for modes 2 and 3 must 
not overlap. Not only will an overlap lead to frequency 
shifts (the jump in amplitude of mode 3 at u) — lj C 2 in 
Fig. [3] comes from such a shift), but we have also ob- 
served that it leads to richer dynamics, including Hopf 



bifurcations with limit cycles |26] . The necessary cri- 
terion for non-overlap can be shown [using Eq. ^} to 
give a lower bound e m i n > 2.2[A/"(xm)]~ 2 Q]~ 1 - Fig. 5b 
shows contours of constant values of e m in<2i- There, re- 
gions close to the edges and the center are excluded. Be- 
cause the responsivity TZ±(r, s) — > 2uji s + 0([1 — s] 2 ) as 
s —> 1, the exclusion of the central area is superficial. 
For example, if we want to use the part of the membrane 
with 0.1 < x, y < 0.9, we have approximately the lower 
bound e > 3Q _1 . For a square membrane of 1 jj,m side 
(M m 760 ag), the present scheme is applicable to masses 
larger than AM m i n « 0.76 ag (assuming Q — 3000). 

Numerical simulations. — To test the scheme we 
implemented an automated mass measurment algorithm 
which numerically integrated the system with a ran- 
domly deposited mass on the membrane. The algorithm 
then determined the frequencies £^1,2,3 and w C 2 iC 3 and cal- 
culated e using Eqs. (|9|) . (fTTT) . and (fl2|) . The results are 
shown in Fig.JSJ;. The relative error in e ranges from 0.1% 
to 98% with the larger errors near the edges where e is 
highly sensitive to position. Masses close to the edges 
could be identified by overlapping responses for modes 2 
and 3 in the nonlinear regimes and were discarded. As can 
be seen, the errors in position of the remaining particles 
are typically small. 

Conclusions. — In conclusion, we have proposed a 
scheme to determine both the position and mass of a sin- 
gle particle adsorbed on a vibrating graphene membrane. 
We have shown that by using bimodal excitation and ex- 
ploiting the nonlinear response of the resonator, measure- 
ments can be restricted to a narrow frequency band near 
the fundamental frequency. Considering that the typical 
resonance frequencies of graphene membranes lie in the 
GHz range, this simplification offers significant experimen- 
tal advantages. These measurements provide information 
about the resonance frequencies and the coefficients of the 
nonlinear terms of the dynamic equations (Kerr constants) 
of the high-order modes. In a resonator without special 
symmetries, the mass and position of the adsorbed particle 
can be determined using the resonance frequency shifts of 
three different modes — measured at a narrow frequency 
band near the fundamental frequency. If the resonator is 
square, it is possible to separate the single-particle adsorb- 
tion events by watching out for changes of the resonance 
frequency of the third mode. Using a gate with a proper 
symmetry, it is possible to determine the mass and posi- 
tion of a adsorbed analyte on the membrane by using the 
resonance frequency shifts of modes 1 and 2 and the fre- 
quencies of the lower-edge bistability regions of modes 2 
and 3. 

As an example we have studied a square membrane with 
an area of 1 /1m 2 , eigenfrequency of 2 GHz and quality 
factor of Q « 3000. For this membrane the sensitivity at 
room temperature (minimum detectable mass change) is 
below 1 zeptogram with a practical operating range in the 
attogram region. This can be compared with, e.g., quartz 



p-5 



J. Atalaya et al. 



crystal microbalances that have mass sensitivities in the 
nanogram range. 
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Appendix. — We here present, for completeness, a 
brief derivation of the location of the bifurcation point 
on the so called backbone curve for the Duffing oscilla- 
tor. Similar derivations can be found in most books on 
nonlinear systems (see for instance |28j). 

Consider a harmonically driven Duffing oscillator x + 
2jx + loqX + kx 3 — po cos(u;i) and introduce slowly in time 
varying action- angle variables r(t) and <f>(t) such that x = 
rsin(wi + cf>) and x — rujcos(ujt + cf>). Substituting these 
expressions into the differential equation and averaging 
over the fast oscillations (see for instance [28]) gives the 
system 



no 



ru>4> 



Po 



uj 2 - up- + (3/i/4)r 2 



Pa 

— cos < 
2 



The frequency response curve is found by solving for the 
stationary regime f = <j> = 0. This amounts to solving the 
frequency response equation 



Aj 2 r 2 cu 2 + r 1 



uj 2 ) + -K,r 2 
4 



Pi 



(13) 



We seek the solution when the bifurcation occur. This is 
exactly the point where ^ = 0. Using this equality while 
taking the derivative with respect to r in the frequency 
response equation (|13[) . leads to an equation for the crit- 
ical frequency u> c (considering here the limit 7 — > 0) for 
transition from the low to large amplitude solution 
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